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ABSTRACT 






The objective of this study was to model and simulate progressive failure 
initiating from a notch tip in a laminated fibrous composite subjected to tensile in- 
plane loading. The micro/macro-level approach was used for this study. The micro- 
level analysis used the 3-D unit cell model while macro-level analysis used the finite 
element analysis technique. A cross-ply laminate with double-edge notches was 
studied to investigate delamination, fiber splitting, and transverse matrix cracking in 
the structure. Numerical results were compared to previous experimental work. 
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I. 



INTRODUCTION 



A. BACKGROUND 

Throughout history technological development has been closely related to our 
ability to use and manipulate materials to suit our needs. Many modem composite 
materials offer design engineers desirable combinations of material properties 
unavailable in more traditional materials. Composites are often used as structural 
members because of their improved strength-to-weight ratio, stiffhess-to-weight ratio, 
toughness, or corrosion resistance. They can also be engineered to provide a particular 
combination of properties by varying such factors as constituent materials, volume 
fraction of constituents, and fabrication/manufacturing processes. The primary drawback 
from most composites is their higher costs. However, increases in component costs are 
frequently offset by decreases in systemic costs due to weight savings. 

B. PREVIOUS WORK 

A number of investigations have been conducted concerning damage and failure 
in composite materials. In the first of their four paper series on Damage Mechanics of 
Composite Materials, Kortschot and Beaumont conducted an experimental study to 
establish a qualitative relationship between the notched strength and terminal damage 
state of double-edge-notched (DEN) cross-ply specimens. They also identified three 
primary failure modes for these specimens: fiber splits in the 0° plies, transverse ply 
matrix cracks in the 90° plies, and triangular delamination zones at the 0/90 interfaces. 
These modes are shown schematically in Figure 1 . They also found that the angles at the 
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apex of the delamination zones measured between 5° and 10° for most specimens. [Ref. 

1 ] 




90° 



Figure 1 Failure Modes in a (90/0)s Laminate— After [Ref. 1] 



The finite element method has been used by many investigators to model and 
predict the response of composite materials to various loads. Chen et al [Ref. 2] and 
Hitchings et al [Ref. 3] both approach the problem using fracture mechanics assumptions 
and strain energy release rate criteria for delamination propagation. Davidson et al [Ref. 
4] simplified that approach further by employing a special crack-tip element in the 
expected damage region. Reedy et al [Ref. 5] modeled delamination by using an eight- 
noded hex constraint element to connect shell elements of different laminae. If their 
failure criterion was met, the connection was considered broken and delamination 
considered to have begun or grown. 

Brewer and Lagace [Ref. 6] approached delamination from the viewpoint of the 
design engineer using known composite materials in structural applications. As such, 
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they sought a quick and efficient method for predicting delamination. They proposed and 
verified a Quadratic Delamination Criterion based on the interlaminar stresses calculated 
efficiently in Kassapoglou and Lagace’s previous work [Ref 7]. The primary advantage 
of their methods is computational simplicity; it enables the designer to predict 
delamination initiation without performing a complete finite element analysis of the 
structure. 

Another approach to failure prediction using the finite element method involves 
using a micromechanical model to separately examine fiber and matrix stresses in 
composite laminae. Ardic et al [Ref. 8] conducted a study examining composites of dis- 
similar laminae using this methodology. Zhu et al [Ref. 9] developed their Direct 
Micromechanics Method (DMM) and compared it with phenomenological results for 
unidirectional composites. 

C. OBJECTIVE 

The objective of this study was to reproduce numerically damage patterns in DEN 
cross-ply continuous fiber composites that had been observed experimentally but which 
are still unexplained analytically. To achieve this a micro-/macromechanical model was 
applied in conjunction with appropriate failure criterion. A parametric study was then 
conducted to investigate the effects of notch size and specimen width on the delamination 
zone size. 
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II. SOLUTION METHODOLOGY 



A. ANALYTIC 

1. Determination of Material Properties 

In order to effectively use composite materials, designers must be able to reliably 
determine their properties including their degradation as a result of damage. Structural 
analysis of composite components is conducted using one of two approaches: micro- or 
macromechanical. In macromechanical analysis, the component is treated as a pseudo- 
homogeneous body. The effective material properties of the pseudo-homogeneous body 
are either predicted through “smearing” or averaging the properties of its constituent 
materials, or they are determined experimentally. The smearing calculation is generally 
based on micromechanical analysis, but after the effective properties have been 
determined, there is no further need for the properties of the constituent materials. 

Micromechanical analysis uses the properties of the constituent materials and 
their interaction to predict the response of the composite material. Three general 
approaches to micromechanical analysis are: the Rule of Mixtures (ROM), semi- 

empirical approaches, and methods of cells. Certain assumptions are common to each of 
these methods; for continuous fiber composites they include: homogeneous, linearly 
elastic, isotropic matrix material; perfect bonding between fiber and matrix; and 
regularly spaced, perfectly aligned, homogeneous fibers [Ref 10]. 

The Rule of Mixtures is a rough tool for predicting composite properties through 
volume averaging; it is based on the strength of materials approach. Longitudinal 
properties result from the assumption that fibers and surrounding matrix material 
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experience the same strain. Transverse properties follow from the assumption that fiber 
and matrix are subjected to the same stress. Experimental work has shown that actual 
longitudinal properties of vmiaxial laminae are comparable to those calculated through 
ROM; calculated transverse properties have proven to be less comparable. Thus, ROM is 
useful as a rule of thumb for estimating properties along the fiber orientation. 

The semi-empirical approach is characterized by the Halpin-Tsai Equations which 
relate composite properties, p, to matrix properties, pm, as follows; 

P 0 ) 

p. 



where the function r\ is of the form shown: 



El 



-1 



7 = 




( 2 ) 



and ^ is determined experimentally and depends on a variety of factors, including fiber 
volume fraction, fiber geometry, loading conditions, and which property, p, is sought. Vf 
is the volume fraction of the fiber. [Ref. 1 1] 

The Rule of Mixtures and the Halpin-Tsai Equations provide means of 
determining smeared properties of a composite lamina in the linear elastic range. To 
progress beyond the linear elastic range into damage analysis and degradation of 
materials as a result of damage, another approach is required. The micromechanical 
method of cells was introduced by Aboudi [Ref. 12]; in it he assumed that a composite is 
a periodic array of matrix and reinforcement (fiber) materials. The analysis then 
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concentrates on a representative unit cell with displacement and traction continuity 
boundary conditions imposed at the interfaces. The unit cell is composed of four 
subcells, one consisting of reinforcement and three consisting of matrix material. 
Although it appears to imply that fibers have square cross-sections, because interfacial 
conditions are imposed on average vice point-wise bases, the geometry is not constrained. 
Pecknold [Ref 13] noted that Aboudi’s model forms the basis for a finite element model 
and conducted an investigation of a simplified unit cell model. Kwon and his colleagues 
have refined the unit cell model to examine stresses in both fiber and matrix components 
at the micromechanical level [Ref 14-20]. This model is particularly suitable for 
analyzing damage, for failure criteria can be applied to both fiber and matrix materials at 
the micromechanical level. This [Ref. 20], then, is the basis for this current w'ork. 

2. Micro-Macro Interaction 

In this work, a combination of micro- and macro-mechanical approaches is used. 
Micro-mechanical analysis is used to determine smeared, effective composite properties. 
Finite element analysis uses the smeared properties in the constitutive equations to 
calculate macro- or elemental stresses and strains. These stresses and strains can then be 
used in a structural analysis or, as in this work, decomposed into micro- or subcell 
stresses and strains. Each subcell consists solely of fiber or matrix material, so failure 
criteria can then be applied to the microstresses and microstrains to determine failure 
initiation, progression, or mode. If failure is found to have occurred, the material 
properties of the constituent materials can then be degraded to reflect that and through the 
smearing process the degradation is carried through to the next iteration of the analysis. 
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Figure 2 illustrates the transitions and exchange of information between the two types of 
analysis. 




Micro/Macro Approach 



STRUCTURE 






Assessment of Failure, Residual strength & stiffness 



Figure 2 Micro/macromechanical Model 

3. The Micromechanical Cell Model 

Although micromechanical models for fibrous composites usually use only four 
subcells, eight subcells were used in this work because the computer code used was 
adapted from one developed for particulate composites which generally use an eight 
subcell model. The principles are the same, but the fiber occupies two subcells along its 
longitudinal direction. The unit cell is a three dimensional solid, rectangular, 
parallelepiped broken down into eight subcells, two containing fiber material and six 
containing matrix material. Taking advantage of symmetry only one quarter of the total 
cell is used; as shown in Figure 3. 
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Figure 3 The Micromechanical Unit Cell 

The size of each subcell is dependent on fiber volume fraction. The stresses and 
strains for each unit cell or element are the volume averaged subcell stresses and strains 
as shown in Equations (3) and (4) 




where cnj and so are the elemental stresses and strains, cr^ and Sy are subcell (k=l,8) 
stresses and strains, and Vf is the fiber volume fraction. Subcells 1 and 2 contain the fiber 
and subcells 3-8 contain matrix material. Stress continuity conditions at subcell 
interfaces are shown in Equation (5). 
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1 2 3 4 5 6 

^11 “”^11 5 ^11 “^11 ? ^11 “^11 

^ ^ _ ^4 ^5 _ _7 

^22 “ ^22 ’ ^22 ^22 ? ^22 “ ^22 

^^^5 2 6 3 7 

^33 ~ ^33 5 ^33 “ ^33 ’ ^33 “ ^33 

^12 “ ^2 5 ^2 “ ^2 5 ^2 

5 6 5 7 6 

^12 ~ ^12 5 ^12 “ ^12 5 ^12 

^ ^ 1—^5 ^5 

^3 “ ^23 5 ^23 “ ^23 5 ^23 

^ _ 4 ^2 _ 6 6 

^23 “ ^23 5 ^23 “ ^23 5 ^23 

^13 “ ^13 5 ^13 “ ^13 5 ^3 

^13 “ ^13 ? ^3 “ ^13 5 <^13 




( 5 ) 



Individual subcells may experience different strains, but the sums of parallel 
longitudinal strains must be equal. A similar compatibility condition is imposed for 
transverse strains. Mathematically these conditions are shown in Equation (6): 




Similar expressions are used for shear strain compatibility. 

The constitutive equations for each subcell follow the form of the generalized 
Hooke’s Law shovm in Equation (7): 

<y^=E“^sl ij,k,l = l,3and« = l,8 (7) 
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The assumption of perfect matrix/fiber bonding can be relaxed if the interface 
between fiber and matrix subcells is modeled using an equivalent spring as shown in 



Figure 4. This changes the strain compatibility equations to those of Equation (8). 
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Figure 4 Micromechanical Unit Cell with Spring Interface 

^ 1 ■•■^1 ~^\ I 1 1 1 



( 8 ) 



Simultaneous solution of Equations (3) through (7) produces smeared composite 
material properties in terms of fiber and matrix material properties and their relative 
compositions. Combining this model with the Finite Element Method detailed in Section 
B, the following sequence of quantities can be calculated: elemental displacements — > 
elemental strains — subcell strains — > subcell stresses — > elemental stresses. 
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In the study of progressive damage, therefore, failure criteria applied to fiber or 
matrix material at the micro-level will result in corresponding degradation of material 
properties at that level. This degradation will then be carried through subsequent 
iterations to reveal the failure processes of the composite. 

B. NUMERICAL 

1. Finite Element Method 

The Finite Element Method of numerical analysis is used to solve the three 
dimensional elasticity problem: to determine stress and strain states of both fiber and 
matrix, and to apply failure criteria for each element. The finite element method converts 
a system of partial differential equations into a larger system of algebraic equations. It 
initially solves for nodal displacement and further processing calculates elemental strains. 
Application of the micromechanical model is used to determine subcell strains and 
stresses; these then serve as entering arguments for evaluating failure criteria for each 
element. 

The derivation of the equations of equilibrium for a three dimensional solid in the 
absence of body forces can be found in any solid mechanics textbook [Ref 21]. Those 
equations are: 
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At first glance there appear to be nine independent variables in these three 
equations; however, application of the moment equilibrium condition reduces this 
number to six: 
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Below is the unit solid element illustrating the sign convention used with stress 
variables. 




Figure 5 Sign Conventions 
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The variables ‘u’, ‘v’, and ‘w’ are used to represent displacements in the ‘x’, ‘y’, 
and ‘z’ directions, respectively. Assuming small displacements, strain is then defined as: 
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These strain-displacement relationships can also be expressed in matrix form as 
shown in Equation (12): 
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The generalized Hooke’s Law relating stress to strain can be expressed as: 

{<t} = [D]{€} (13) 



where {a} is the 6x1 stress column vector, {s} is the 6x1 strain vector and [D] is a 6x6 
matrix of material properties. Thus Equations ( 9)-(13) combine as fifteen equations with 
fifteen unknowns, including several partial differential equations. 

The derivation of the finite element equations from the three dimensional 
elasticity equations above is based on course notes and commonly available current 
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textbooks [Ref. 22]. The method of weighted residuals as modified by Galerkin is 
employed to develop the displacement based finite element equations. Eight noded 
isoparametric rectangular parallelepiped elements are used for the formulation. Linear 
shape functions are used for analytical and computational simplicity. Detailed derivation 
according to these precepts is continued below. 

The first step in converting the partial differential equations into algebraic 
equations is to apply the method of weighted residuals to the three dimensional stress 
equilibrium Equations (9). To this end, each of the three equations of Equation (9) are 
multiplied by a weight function, Wj (i=l, 2, 3), which is continuous over the physical 
domain of the problem. Then, the three new product equations are integrated over the 
entire problem domain. The goal is to chose weight functions Wj which are orthogonal to 
the initial residuals of the equilibrium equations such that the integral of their product is 
zero. If ‘V’ is the domain volume of the problem, the weighted residual equations are 
shown in Equation (14): 



At this point it should be noted that boundary conditions must be specified for Equation 
(9). The boundary conditions may be specified as either (a) essential or geometric 
boundary conditions in which some surface displacements are specified, (b) natural or 




(14) 
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stress boundary conditions in which surface tractions are specified, or (c) a combination 
of these types of boundary conditions. Before applying boundary conditions, further 
manipulation of the weighted residuals is required as shown by Equation (15). 
Integration by parts will be performed in order to weaken the continuity requirements on 
the approximate solution. 



A = 111 - 

V _ 

= III’- 

V , 

^3 = in'- 



(J h T h T 



^2 SV, 

sv, M, m. ' 



dV + + r,MW,)lS 

s 
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dV + 



In Equation (15) ‘S’ is the domain surface boundary, ‘V’ is the domain volume, and ‘nx’, 
‘ny’, and ‘nz’ are outward unit normal directional cosines in the ‘x’, ‘y’, and ‘z’ 
directions, respectively. The boundary stress conditions are defined by surface tractions 
as shown in Equation (16): 

<Px = + ^xA + ^.z«- 

S=Tn+an+Tn 

Ty "xy'‘x '-'y'y "yz'‘z 

+ Ty,ny + a.n. 

Equations (15) and (16) can be combined and expressed in matrix form as: 
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Equation (17) can be further modified by separating the column vector inside the volume 
integral into the product of a 3x6 matrix and a 6x1 vector. This step, shown in 
Equation(lS), isolates the weight function derivatives in the matrix from the stress terms 
in the vector. 
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Substitution of the strain-displacement relation of Equation (11) into the 
constitutive relations given by Equation (13) provides another useful identity as shown in 
Equation (19). 
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Equation (12) is now substituted for the column vector on the left hand side of Equation 
(19). Equation (18) and the modified Equation (19) are substituted into Equation (17) to 
produce Equation (20). 
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( 20 ) 



The preceding manipulations of the elasticity relations have isolated displacement 
as the desired quantity. Now, to discretize the problem for an algebraic computational 
solution, assume that over a small domain each of the displacements can be represented 
by a polynomial. Discretization divides a three dimensional domain into small-but-finite 
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volume elements. In this formulation each finite element has eight nodes; each node can 
have three orthogonal displacements. Each element, therefore has a total of twenty-four 
degrees of freedom (dof). Calling ‘Uj’, ‘vi’, and ‘Wj’ the displacements at each node, 
elemental displacements can be defined in terms of polynomial shape functions as 
follows: 

u = Y,H,u,; v = w = Y,H,w. 

/=! ;=I i=I 



The first partial derivatives of the shape functions exist as shown in Equations (22). 
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The 3x1 elemental displacement vector can be expressed as the product of a 3x24 shape 
function matrix and a 24x1 nodal displacement vector. The product of the 6x3 partial 
differential matrix of Equation (20) and the 3x24 shape function matrix are typically 
combined into one 6x24 matrix. This matrix is referred to as the ‘B’ matrix in this 
formulation. In the shorthand notation shown in Equation (23), the ‘B’ matrix can be 
partitioned into sub-matrices, ‘Bj’, where i=l to 8. 

W=P.][^][^][Sd[Ss]M[S>][All <23) 

In this notation, the sub-matrices are defined in Equation (24). 
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where i = 1 to 8 



The Galerkin method takes the weighting functions as equivalent to the shape 
functions as shown in Equation (25). 
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The Galerkin method only requires that the weighting function be continuous over small 
discrete intervals which correspond to the sides of the finite elements. Based on the 
Galerkin weighting functions shown above, the weighting function matrix of Equation 
(20) can now be written in terms of the shape functions: 
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Equation (26) can be incorporated in Equation (20) as follows: 
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In Equation (27), {d} is the displacement vector which has grown from the 3x1 
vector of ‘u’, ‘v’, and ‘w’ to a 24x1 vector of the nodal displacements ‘uj’, ‘vj’, and ‘Wj’. 
The 3x1 vector in the surface integral on the right-hand side has an 8x1 vector in each 
term so that the entire integrand is shorthand for a 24x1 vector of discretized surface 
boundary tractions. The original three equilibrium elasticity partial differential equations 
have now been transformed to a matrix equation in 24 terms for each solid element. The 
integrals in Equation (27) may be easily solved if simple shape functions are chosen and 
if the modeled geometry is relatively simple. 

If the modeled geometry is complex, the finite element geometry can be 
simplified by use of a transformation mapping to another reference space. In order to 
map ‘x’. ‘y’, and ‘z’ coordinates of an irregularly shaped element onto ‘r’, ‘s’, and ‘t’ 
coordinates for a rectangular parallelepiped, a Jacobian transformation matrix is required. 
The Jacobian transform matrix scales each component of the original space to a new 
space. Equation (28) shows the Jacobian for an ‘xyz’ system mapped to an ‘rst' system: 
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In terms of the shape functions and nodal points the Jacobian is shown in Equation (29). 
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Since the finite element integral equation includes the partial derivatives of the shape 
functions with respect to the ‘xyz’ coordinate system, the inverse of the Jacobian is 
required. Let [J]’'=[r], where [F] is a 3x3 matrix. The shape function derivatives with 
respect to the ‘xyz’ system are now expressed in Equation (30) with respect to the ‘rst’ 
system. 
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X 

Equation (30) is used to calculate the [B] and [B] matrices in the ‘rst’ system. Equation 
(29) is used to transform the volume differential. Application of these transformations to 
the left-hand side of Equation (27) is shown in Equation (31). 
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( 31 ) 



V y 

The transformation to the ‘rst’ coordinate system results in simplified finite 
element integrals. The resultant transformed elements are termed isoparametric elements. 
The term isoparametric refers to the equality between the degree of the equations for 
transforming ‘xyz’ into ‘rst’ coordinates and the degree of the shape functions for 
estimating displacement. The shape functions in the ‘rst’ system are shown in Equation 
(32). 

/X=-50+'-)0-jXi+') 

H,=4(i+'-Xi-s)0-<) 

H,=-l(l+r)0+sXl-0 

//,=i(i-/-Xi+^)0+/) 

fli=-i0-''Xi-^)0+') 

ff,=-S(i-'-Xi-j)(i-') 

Gauss Quadrature is used to numerically calculate the integrals. The volume 
integral is redefined as the triple summation from 1 to the number of integration points 
(NIP) of the integrand evaluated at Gauss integration points (r„ Si, and t;) multiplied by 
weight factors as shown in Equation (33). 
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NIP NIP NIP 



(33) 



The results of the numerical integration may vary over the elements in the domain of the 
model. Each of these results can be expressed as a 24x24 elemental stiffness matrix, 



Recall the surface traction boundary conditions from the right-hand side of Equation (20). 
The integration of the directional component of the applied stress over the element 
surface area to which the stress is applied is equal to the external force applied to the 
element. The result of the integration is a 24x1 vector {Fe} which is equivalent to the 
force in the ‘x’, ‘y’, and ‘z’ direction applied to a surface of the solid. Substituting these 
resultant terms into Equation (20) completes the finite element derivation shown in 
Equation (34). 



The finite element method involves combining these elemental matrix equations 
with given boimdary conditions to form a large system of simultaneous equations to be 
solved numerically. 



The actual finite element analysis for this study was conducted using a 
FORTRAN driver for three-dimensional static analysis based in large part on the 
subroutines developed by Akin [Ref. 23]. Pre- and post-processing were conducted using 
various MATLAB programs. 



[K.]. 






(34) 



2. Implementation 
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The preprocessor generated an appropriately formatted text file containing 
material properties, nodal coordinates, element connectivity, boundary conditions, and 
loading data. 

The FORTRAN program used several subroutines to read the input data and store 
it semi-dynamically. All data was stored in one of two massive arrays, one for integers 
and one for floating point variables, and pseudo-pointers were used to track the location 
of the first element of each matrix of interest. The memory allocation for these two 
arrays was calculated based upon the input data. After the input and housekeeping 
subroutines were completed, the heart of the analysis was controlled by subroutine 
FEMST. It called subroutine EL3D8A to calculate elemental stiffness matrices and 
assemble them into the system stiffness matrix. The stiffness matrices were calculated 
using the micromechanical model and Gauss Quadrature with two integration points in 
each coordinate direction for a total of eight integration points per element. Subroutine 
MA3D6F calculated material property matrices for composite elements, including 
transforming those matrices to account for fiber orientation. After the system matrix was 
assembled, boundary conditions were applied, and Equation (34) was solved for the nodal 
displacements using subroutines FACTOR and SOLVE. FACTOR uses the Cholesky 
method to factor the square system stiffness matrix into the product of a lower triangular 
matrix and its transpose. SOLVE uses Cholesky-Gauss methods of forward and back 
substitutions to complete the solution. 

In the process of calculating the residual, elemental stresses and strains are 
calculated and then decomposed into subcell stresses and strains. These values can then 
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be used to evaluate matrix and fiber failure criterion and degrade the appropriate 
corresponding material properties if necessary. 

The FORTRAN program writes its output to five main files: OUTPUT, DISPL, 
MATSTS, FIRSTS, and STRESS. OUTPUT is essentially an echo of the input file; 
DISPL is the u, v, and w displacements of each node; MATSTS contains the matrix 
subcell stresses at each integration point; FIRSTS contains the fiber subcell stresses at 
each integration point; and STRESS contains the complete stress state of each element. 
These data files were then read into various MATLAE programs for further analysis and 
visualization. 

C. FAILURE CRITERION 

The nature of the failure modes being investigated in this study, delamination and 
splitting, dictates that matrix material behavior be of primary interest. These phenomena 
are of interest because they can cause a structure to fail at much lower than expected 
loads. Delamination is an inter-laminar process; as such, the inter-laminar, through- 
thickness, or z-direction stresses are expected to be the controlling factor in its initiation. 
Rrewer and Lagace [Ref. 6] proposed the following criterion for delamination: 
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Where 7} is the shear strength, 7} and 7^ are normal strengths in tension and 
compression, respectively; and Txz and Xyz are interlaminar shear stresses, a* is 
interlaminar tensile stress, and a' is interlaminar compressive stress. This criterion will 
be the starting point of this investigation in which the failure criterion will be applied to 
generate a delamination zone similar to that observed experimentally. Assuming that the 
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isotropic matrix material is twice as strong in tension as it is in shear and four times as 
strong in compression as it is in tension, Equation (35) can be rearranged in terms of a 
single strength value as shown in Equation (36); 




(36) 



In this form, this criterion can be evaluated without necessarily knowing the 
strength of the material in question. A stress profile calculated from the left-hand side of 
Equation (36) can be plotted and used to evaluate the suitability of the criterion. 

Similar criterion was applied using stresses transverse to the fiber direction to 
evaluate fiber splitting and transverse cracking. 
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III. RESULTS AND DISCUSSION 



A. COMPARISON OF 2-LAYER MODEL TO 3-LAYER MODEL TO 
EXPERIMENTAL RESULTS 

In their experimental work Kortschot and Beaumont [Ref. 1] applied tensile 
loading to double-edge-notched graphite-epoxy specimens of various cross-ply lay-ups 
and sizes. They observed that delamination along the 0/90 interfaces progressed from the 
notch tip in a triangular zone upward and toward the center of the specimen. As the splits 
in 0° plies grew, the aspect ratio of the delamination triangle remained constant. The 
measured angles at the apex of each specimen’s delamination triangle were found to 
consistently remain between 5° and 10°. 

In this study two models were used to simulate the behavior observed by 
Kortschot and Beaumont. A two-layer model consisting of two laminae, one with fibers 
oriented along the direction of loading (0°) and one with fibers oriented transverse to the 
loading direction (90°), was compared with a three-layer model consisting of the same 
two laminae plus a thin “delamination layer” consisting solely of isotropic matrix 
material inserted between the two laminae. The same material properties and geometry 
were used for each model. The material properties of the tu'o laminae of each model 
differed only by their orientation. The properties of the constituents of the specific 
composite used in the experimental work were unavailable, but values were chosen so 
that they were within approximate ranges for graphite fibers and epoxy matrix and that 
their smeared values were comparable to those cited by those investigators; the values 
used are shown in Table 1. The lay-up modeled was (90/0)s, that is, two 0° plies 
sandwiched between a pair of 90° plies. Geometric and laminar symmetry conditions 
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allowed for a proper model of the specimens using only one-eighth of the original. 
Planes of symmetry were the mid-plane both longitudinally and transversely and the 
plane between the 0° plies. Nodes located on these planes were constrained in the 
direction(s) normal to their respective planes of symmetry. As in the experimental work, 
the modeled notch angle was 60°. 



Table 1 Lamina Dimensions and Material Properties 



Width 


5mm 


Length 


10mm 


Thickness 


1.5mm 


Notch Length/Specimen Width (2a/w) 


0.125 


Fiber Longitudinal Young’s Modulus (Efi) 


202.5 GPa 


Matrix Young’s Modulus (Em) 


3.45 GPa 


Fiber Transverse Young’s Modulus (Eft) 


15.75 GPa 


Fiber Poisson Ratio (12 direction) 


0.29 


Fiber Poisson Ratio (23 direction) 


0.25 


Matrix Poisson Ratio 


0.35 


Fiber Shear Modulus (12 direction) 


15.75 GPa 


Volume Fraction 


0.66 



To evaluate the suitability of the two-layer model, the left-hand side of Equation 
(36) was applied to the volume averaged matrix subcell stresses for each element in the 
0° lamina. This profile is shown for the entire modeled section in Figure 6 and for the 
refined mesh area in Figure 7. 
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Figure 6 Delamination Zone 0° Ply 2 Layer Model 
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X 1 0 Quadratic Criterion 0 degree Layer 2 Layer Model 
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Figure 7 Delamination Zone 0° Ply 2 Layer Model 
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The same sort of profile was generated for the three layer model based on the 
stress states of the elements in the delamination layer. These are shown in Figure 8 and 
Figure 9. 




Figure 8 Delamination Zone 3 Layer Model 
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Figure 9 Delamination Zone 3 Layer Model 



It is apparent from Figures 6-9 that the three-layer model more closely resembles 
the delamination zone observed experimentally. The apex angle of the delamination 
triangle (superimposed in Figure 9) in the three layer model is approximately 10°, 
comparing favorably with experimental observations. While the two layer model’s stress 
profile does not have the characteristic triangular shape of a delamination region, a 
similar profile generated from its transverse matrix subcell stresses does compare 
favorably with observed 0° ply splitting. This profile is shown in Figures 10 and 1 1. 
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Figure 10 Fiber Splitting in 0° Ply, 2 Layer Model 
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Figure 1 1 Fiber Splitting in 0° Ply, 2 Layer Model 



The three layer model not only displays a more realistic delamination zone, the 
transverse stress profiles in the 0° and 90° laminae also reflect fiber splitting and 
transverse cracking, respectively. As observed experimentally and expected with 
graphite fibers, transverse cracking permeates the 90° lamina. These effects are 
illustrated in Figures 12-14. 
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Figure 12 Fiber Splitting in 0° Ply, 3 Layer Model 
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Figure 13 Fiber Splitting in 0° Ply, 3 Layer Model 
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Figure 14 Transverse Cracking 90° Ply, 3 Layer Model 
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PARAMETRIC STUDY 



A parametric study varying notch length and specimen width was conducted to 
further evaluate the quadratic delamination criterion as applied to the three layer model. 
Using the original crack length to width ratio (2a/w=0.125), samples 5, 10, 20, and 40 
mm wide were simulated. The resulting delamination zones are shown in Figures 15-18. 
The apex angles of those zones ranged from 6°-12°. 

In the pre-processor used, the nodes were placed as a function of crack length and 
specimen width and length; as such, each of these different specimens has the same 
number of elements. Therefore, the wider specimens have larger elements and some 
detail is lost. Nevertheless, the influence of the fiber splits and the general triangular 
shape of the delamination zones can be seen. 



Delamination Zone 2a/w=0.125 
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Figure 15 Delamination Zone w=5mm 
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Figure 16 Delamination Zone w=10mm 
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Figure 1 7 Delamination Zone \v=20mm 



41 




Delamination Zone 2a/w=0.125 
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Figure 1 8 Delamination Zone w=40mm 



For the same four widths crack length to width ratio was also varied. Figures 19- 
22 show delamination zones for those specimens. For the first three of these specimens 
apex angles range from 6°-13°. Once again the mesh effects appear to suggest that 
delamination is a size-dependent phenomenon. This coarsening is most apparent in the 
40mm specimen. Figure 22. 
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Delaminaton Zone 2aAv=0.25 
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Figure 19 Delamination Zone w=5mm 
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Figure 20 Delamination Zone w= 10mm 
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Delamination Zone 2a/w=0.25 
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Figure 21 Delamination Zone w=20mm 



Delamination Zone 2a/w=0,25 




Figure 22 Delamination Zone w=40mm 
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IV. CONCLUSIONS AND RECOMMENDATIONS 



The combined micro- and macromechanical approach to analysis of composite 
materials is an effective method of exploring and simulating the initiation of damage. It 
efficiently enables the incorporation of constituent homogeneous material behavior into 
analysis of heterogeneous materials. In this work micro-/macromechanical analysis was 
used to simulate each of the three experimentally observed modes of failure for DEN 
cross-ply fibrous composites. A thin layer of homogeneous matrix elements was 
necessary to properly model observed delamination behavior. The delamination layer 
made the transmission of stresses and damage across the laminar boundaries conveniently 
observable. The quadratic delamination criterion applied to delamination layer elements 
did correctly predict delamination zone shape for a variety of geometries. Proper mesh 
refinement is necessary to correctly simulate specimen behavior in the vicinity of the 
notch tip. 

Future research in this area should include both experimental and simulated 
examinations of multi-axial loading. This work examined failure modes that occurred 
exclusively in the matrix regions; incorporation of fiber failure and degradation should be 
included in future work. 
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